Os índigenas Pima - diabetes

Gabriel de Jesus Pereira

2025-02-19

Sobre o banco de dados

Sobre o banco de dados

  • Foi originado do National Institute of Diabetes and Digestive and Kidney Diseases.

  • Essa base foi criada com base nos Pima, um grupo de nativos americanos que vivem em uma área que atualmente abrange o centro e o sul do estado do Arizona.

  • O objetivo é prever se um paciente tem diabetes com base em algumas características do paciente.

Descrição das colunas

  • Pregnancies: Número de gestações da paciente.

  • Glucose: Concentração de glicose plasmática em teste oral de tolerância à glicose.

  • BloodPressure: Pressão arterial diastólica (mm Hg).

  • SkinThickness: Espessura da dobra cutânea do tríceps (mm).

  • Insulin: Nível sérico de insulina (mu U/ml).

  • BMI: Índice de Massa Corporal (peso em kg / altura em m², IMC).

  • DiabetesPedigreeFunction: Histórico familiar de diabetes (mede a predisposição genética).

  • Age: Idade do paciente.

  • Outcome: 0 é não diabético e 1 é diabético.

Análise exploratória de dados

Quantidade de valores ausentes

library(kknn)
library(naivebayes)
library(discrim)
library(tidyverse)
library(tidymodels)
library(visdat)
library(corrplot)

data <- read_csv("diabetes.csv") |> mutate(Outcome = as.factor(Outcome))
cols <- c("Glucose", "BloodPressure", "SkinThickness", "Insulin", "BMI")
data[cols][data[cols] == 0] <- NA


Ausentes <- data |>
  vis_dat() +
  labs(
    y = "Observações"
    ) +
  scale_fill_manual(
    values = viridis::mako(10)[5:6],
    labels = c("Qualitativa", "Quantitativa", "Ausentes")
    ) +
  theme(
    legend.position = "bottom",
    legend.title = element_blank(),
  )

Quantidade de valores ausentes

Apenas as colunas de glicose, pressão arterial, espessura da dobra cutânea do tríceps, nível de insulina e IMC. Aquela que mais possui valores ausentes é a variável de nível de insulina.

Frequência

set.seed(123)
data_split <- initial_split(data, prop = 0.8, strata = Outcome)
train_data <- training(data_split)
test_data <- testing(data_split)

rec <- recipe(Outcome ~ ., data = train_data) %>%
  step_impute_median(all_numeric_predictors())

prep_rec <- prep(rec)
train_data_eda <- bake(prep_rec, new_data = train_data)
test_data_eda <- bake(prep_rec, new_data = test_data)

Frequência <- train_data_eda |>
  ggplot(aes(x = forcats::fct_infreq(as.factor(Outcome)))) +
  geom_bar(
        aes(
          y = after_stat(count) / sum(after_stat(count)),
          fill = Outcome
          ),
        show.legend = FALSE,
  ) +
  theme_bw() +
  theme(axis.title.x = element_blank()) +
  scale_fill_manual(
    values = viridis::mako(5)[3:5]
    ) +
  labs(y = NULL, x = NULL)

Frequência

A base de dados possui um problema de desbalanceamento das classes. A classe com maior frequência é de não diabéticos, com pouco mais de 60%.

Distribuição

Distribuição <- train_data_eda |>
  pivot_longer(
    !Outcome,
    names_to = "variable",
    values_to = "values"
  ) |>
  mutate(Outcome = as.factor(Outcome)) |>
  ggplot(
    aes(
      x = values,
      fill = Outcome
      )
    ) +
  geom_violin(
    aes(y = variable),
    draw_quantiles = c(0.25, 0.5, 0.75),
    ) +
  labs(x = "", y = "") +
  facet_wrap(vars(variable), scales = "free") +
  scale_fill_manual(
    values = viridis::mako(5)[3:5],
    name = "Outcome"
    ) +
  theme_bw() +
  theme(axis.text.y = element_blank())

Distribuição

As variáveis apresentam uma assimetria negativa, sendo a variável de pressão arterial a única mais simétrica. Além disso, observando a distribuição dos diabéticos e não diabéticos para cada variável, eles apresentam um padrão de comportamento bastante similar.

Correlação

data_cor <- train_data_eda |>
  (\(x) {
     scaled <- x |>
       select_if(is.numeric) |>
       scale() |>
       as_tibble()

     x |>
       select_if(is.character) |>
       cbind(scaled)
  })()

data_cor |>
  select_if(is.numeric) |>
  cor() |>
  corrplot(
      method="circle",
      col = viridis::mako(200),
      diag = FALSE,
      addgrid.col = "#00000040",
      tl.pos = "l",
      tl.cex = 0.6,
      tl.col = "black",
      tl.offset = 0.25,
      cl.length = 9,
      cl.cex = 0.6,
      cl.ratio = 0.25,
      family = "mono"
  )

Correlação

A maioria das variáveis independentes não apresentam uma alta correlação. No entanto, os níveis de insulina tem uma alta correlação com a glicose, o IMC tem uma alta correlação com a espessura da dobra cutânea do tríceps e a idade tem uma alta correlação com a quantidade de gravidez.

Preparação dos dados

Preparação dos dados

data_split <- initial_split(data, prop = 0.8, strata = Outcome)
train_data <- training(data_split)
test_data <- testing(data_split)

rec <- recipe(Outcome ~ ., data = train_data) |>
  step_impute_median(all_numeric_predictors()) |>
  step_mutate(across(all_numeric_predictors(), log1p)) |>
  step_mutate(
    N0 = BMI * SkinThickness,
    N8 = Age / Insulin,
    N13 = Glucose / DiabetesPedigreeFunction,
  ) |>
  step_nzv(all_numeric_predictors()) |>
  step_normalize(all_numeric_predictors())

Preparação dos dados

  • Os dados foram divididos entre teste e treino, com 20% para o conjunto de teste e estratificado pela variável dependente.

  • Os valores ausentes foram substituídos pela mediana das variáveis.

  • Foram criadas novas variáveis, idade por nível de insulina, glicose por predisposição de ter diabetes e o produto entre IMC e a espessura da dobra cutânea do tríceps.

  • Todas as variáveis foram transformadas com \(\log \left(1 + x\right)\) e normalizadas.

Algoritmos e validação cruzada

cv_folds <- vfold_cv(train_data,
                     v = 10,
                     strata = Outcome
                     )

knn_spec <- nearest_neighbor(neighbors = tune(), weight_func = tune()) |>
  set_engine("kknn") |>
  set_mode("classification")

nbayes_spec <- naive_Bayes(smoothness = tune(), Laplace = tune()) |>
  set_engine("naivebayes") |>
  set_mode("classification")

lda_spec <- discrim_linear() |>
  set_engine("MASS") |>
  set_mode("classification")

qda_spec <- discrim_quad() |>
  set_engine("MASS") |>
  set_mode("classification")

logistic_spec <- logistic_reg() |>
  set_engine(engine = "glm") |>
  set_mode("classification")

wf = workflow_set(
  preproc = list(rec),
  models = list(
    knn_fit = knn_spec,
    nbayes_fit = nbayes_spec,
    lda_fit = lda_spec,
    qda_fit = qda_spec,
    linear_fit = logistic_spec
    )
  )  |>
  mutate(wflow_id = gsub("(recipe_)", "", wflow_id))

grid_ctrl = control_grid(
  save_pred = TRUE,
  parallel_over = "resamples",
  save_workflow = TRUE
)

grid_results = wf  |>
  workflow_map(
    seed = 123,
    resamples = cv_folds,
    grid = 50,
    control = grid_ctrl
  )
autoplot(grid_results, metric = "roc_auc") +
  theme_bw() +
  labs(y = "Métrica", x = "")
autoplot(grid_results, select_best = TRUE) +
  theme_bw() +
  labs(y = "Métrica", x = "")

Algoritmos e validação cruzada

  • A validação cruzada foi realizada com 10 folds, estratificado pela variável dependente. Assim, foram otimizados os algoritmos a quantidade de vizinhos e a métrica de distância do algoritmo de k-vizinhos mais próximos. Foram otimizados, para o naive bayes, o parâmetro de suavidade para o limite da classe e a correção da suavidade chamado de Laplace.

Resultado geral da otimização

Resultado na base de teste

results_acc = workflowsets::rank_results(grid_results,
                                         select_best = TRUE,
                                         rank_metric = "roc_auc") |>
  filter(.metric == "roc_auc") |>
  dplyr::select(wflow_id, mean, std_err, model, rank)

best_set_linear = grid_results %>%
  extract_workflow_set_result("linear_fit") %>%
  select_best(metric = "roc_auc")

best_set_knn = grid_results %>%
  extract_workflow_set_result("knn_fit") %>%
  select_best(metric = "roc_auc")

best_set_nbayes = grid_results %>%
  extract_workflow_set_result("nbayes_fit") %>%
  select_best(metric = "roc_auc")

best_set_lda = grid_results %>%
  extract_workflow_set_result("lda_fit") %>%
  select_best(metric = "roc_auc")

best_set_qda = grid_results %>%
  extract_workflow_set_result("qda_fit") %>%
  select_best(metric = "roc_auc")

test_results <- function(rc_rslts, fit_obj, par_set, split_obj) {
  res <- rc_rslts %>%
    extract_workflow(fit_obj) %>%
    finalize_workflow(par_set) %>%
    last_fit(split = data_split,
             metrics = metric_set(
              accuracy,roc_auc,
              brier_class))
  res
}

test_results_linear = test_results(grid_results, "linear_fit", best_set_linear,df_split)
test_results_knn = test_results(grid_results, "knn_fit", best_set_knn,df_split)
test_results_nbayes = test_results(grid_results, "nbayes_fit", best_set_nbayes,df_split)
test_results_lda = test_results(grid_results, "lda_fit", best_set_lda,df_split)
test_results_qda = test_results(grid_results, "qda_fit", best_set_qda,df_split)

metrics_table = rbind(
  collect_metrics(test_results_linear)$.estimate,
  collect_metrics(test_results_knn)$.estimate,
  collect_metrics(test_results_nbayes)$.estimate,
  collect_metrics(test_results_lda)$.estimate,
  collect_metrics(test_results_qda)$.estimate
)
metrics_table <- round(metrics_table, 4)
rnms = c("logistic","k_nn","naive_bayes", "lin_discr", "quad_discr")
metrics_table <- cbind(rnms, metrics_table)
metrics_table <- metrics_table %>% as_tibble()

colnames(metrics_table) = c(
  "Modelo", "accuracy", "roc_auc", "brier")
metrics_table |>
  arrange(desc(roc_auc)) |>
  knitr::kable()
Modelo accuracy roc_auc brier
lin_discr 0.7857 0.8652 0.1443
logistic 0.7987 0.8546 0.1484
k_nn 0.7468 0.8326 0.1586
quad_discr 0.7727 0.8238 0.1663
naive_bayes 0.7013 0.8128 0.1842

Resultado na base de teste

O modelo que obteve o melhor resultado foi o discriminante linear, com base na curva ROC e o brier. Na base de teste, ele obteve a curva ROC foi de 0,8652, o brier 0,1443 e obteve uma acurácia de 78,57%. O que pior perfomou foi o naive bayes, obtendo as piores estatísticas para todas as métricas consideradas.